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I. INTRODUCTION 

The Boltzmann equation (BE) lies at the basis of classical 
and quantum kinetic theory of gases. It provides a detailed 
picture of the time evolution of a dilute gas towards a ther- 
mal equilibrium state, which constitutes the essence of the 
H-theorem. This celebrated result gave rise, historically, to 
the first clear insurgence of irreversibility into deterministic 
equations of motion. The nonlinear integro-differential nature 
of the BE prevented, so far, an exact solution. Perturbative 
methods and kinetic toy models have been devised such to 
get partial answers. The Chapman-Enskog expansion (CE) 
was, in particular, the first important success in this direction 
[1], as it allowed to consistently derive hydrodynamics laws 
from their microscopic counterpart and to obtain rigorous ex- 
pressions for transport coefficients. The CE method is based 
upon a perturbative expansion of the distribution function in 
terms of the Knudsen number e, defined as the ratio between 
the mean free-path and a macroscopic hydrodynamic length. 
This is supposed to be a "smallness" parameter, in that the 
series converges only for e — > 0. By increasing the order of 
the expansion one should not expect to capture larger extents 
of the "true" solution of the BE, since, as it was pointed out 
by Bobylev [2], one has to face divergencies of the acous- 
tic modes in the dispersion relation, which are inherently re- 
lated to the procedure of truncation. In order to tackle this un- 
physical feature of post-Navier Stokes hydrodynamics, some 
regularization methods were borrowed from functional anal- 
ysis in order to restore the H-Theorem [3]. Another route, 
which attempts a non-perturbative approach to solve the BE, 
is based upon the notion of Invariant Manifold [4]. Through 
this method, one assumes a priori a separation of the hydro- 
dynamic time scale and the kinetic time scale and postulates 
the existence of a stable Invariant Manifold (IM) in the space 
of distribution functions, which is parameterized with the val- 
ues of the hydrodynamic fields: particle number, velocity, and 
temperature. In this paper we address the study of the spec- 
trum of hydrodynamic excitations in a Maxwell gas, employ- 
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ing the latter non-perturbative approach, which will allow to 
find exact transport coefficients at arbitrary length scales. The 
paper is organized as follows: in Sec. II we review the eigen- 
value problem associated with the linearized Boltzmann equa- 
tion and recall that hydrodynamic modes at finite wavevec- 
tor can be obtained as eigenvalues of a perturbed linear op- 
erator. Next, in Sec. Ill, we motivate and derive the invari- 
ance equations (details collected in Appendix A) and consider 
the case of Maxwell molecules (Sec. Ill A), whose associated 
eigenvalue problem for the unperturbed operator is analyti- 
cally solvable (Appendix B). Postulating the existence of an 
IM, we solve the eigenvalue problem for arbitrary wavevectors 
and find generalized transport coefficients which recover the 
Green-Kubo formulas in Sec. Ill B. Further, in Section IV, we 
determine the spectrum of density fluctuations and formulate 
a hypothesis about the features of finite wavelengths hydrody- 
namics. Conclusions are drawn in Sec. V. 



II. EIGENVALUE PROBLEMS FOR THE BOLTZMANN 
EQUATION AND HYDRODYNAMICS 

The dynamics of the fluctuations of hydrodynamic fields 
(particle number, momentum, temperature) as induced by the 
properties of the underlying microscopic or kinetic equation, 
is an important issue in statistical mechanics which dates back 
to the seminal work by Onsager [5]. In this section we focus 
upon the BE and show in a general setting how it features 
equilibration through some generalized frequencies (inverse 
of characteristic collision times). The way how these general- 
ized frequencies give rise and affect the decay rates of some 
collective fluctuations (hydrodynamic modes) of the macro- 
scopic fields, is still an issue which lacks a rigorous foun- 
dation. The reason is that the hydrodynamic equations, as 
derived from the BE, are not closed and hence, some (semi- 
phenomenological) approximations for higher order moments 
need to be included. In particular, the celebrated Navier- 
Stokes-Fourier (NSF) approximation was the first historically 
relevant attempt in this direction. We start by reviewing some 
results formerly obtained by Resibois [6]. He shed some pre- 
liminary light upon the connection between the generalized 
frequencies and the hydrodynamic modes by solving, via per- 
turbation theory, the eigenvalue problems associated indepen- 
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dently to the BE and to the NSF equations of hydrodynamics. /(k, c, u) — e u ' t f(k 7 c, t)dt, where w defines a com 



A. The Boltzmann equation 



The BE reads: 



dtf = -y-Vf + Q[fj], 



(1) 



where Q denotes a non-linear integral collision operator. We 
introduce the thermal velocity vt = \j2kBTojm, the dimen- 
sionless peculiar velocity c = (v — Uq)/vt and the equi- 
librium values of macroscopic fields: equilibrium particle 
number no, equilibrium mean velocity Uo = 0, and equilib- 
rium temperature To. The global Maxwellian is defined as: 
/ GM = (n /v T )f (c) where f (c) = 7r~ 3 / 2 e _c denotes a 
Gaussian in velocity space (c = |c|). We consider only small 
disturbances from the global equilibrium. After passing over 
to Fourier space, we write the distribution function (cf. also 
Tab. I) as: 



/(M,t) = / 



LM 



Sf, 



(2) 



where / LM denotes the local Maxwellian to be made precise 
in Sec. Ill, and 5f the deviation from local equilibrium. An 
alternative notation is introduced via Sf = f GM Sip. Consider- 
ing a co-moving reference frame and linearizing the collision 
operator around global equilibrium, one obtains from (1) 



Vt 



-d t f 



ik ■ cf + LSf, 



Vt 



-L, 



(3) 



where we made use of the fact that L/ LM = 0. The linearized 
Boltzmann collision operator, L, assumes the form 

LSf = J J dndc 1 a(n,g)gf M (c 1 ) x 

[<V(M) + MMO - <V(M0 - ^(k,ci].(4) 

Here, a(Sl,g) is the scattering cross section, g = \v — Vi\, 
and v, vi are the velocities of the particles entering the bi- 
nary collision. In the remainder ot this section, we will focus 
our attention upon the operator A = L — ik ■ c, whose spec- 
tral properties determine the time evolution of the distribution 
function. This is readily seen by considering the Laplace time 
transform of Eq. 3 (to be further discussed in Sec. IV) and by 
inspection of the inverse transform, which reads as: 



/(k,c,i) 



1 



2m J (z - A) 



<te/(k,c,0), 



(5) 



where the closed path encircles all the poles of the integrand 
function. Through the Spectral Theorem, we regard these 
poles as coinciding with the spectrum of the operator A. The 
investigation of the spectral properties of such an operator is, 
in fact, a longstanding issue in kinetic theory [7]. In order 
to study the eigenvalue problem associated with A, we intro- 
duce the Fourier time transform of the distribution function 



plex valued quantity. Then, (3) reduces to: 

A/ = uf, 



(6) 



which constitutes the starting point of our analysis. 

In the present paper, functions / M = / GM (^ M will be re- 
garded as vectors in a Hilbert space, whose scalar product is 
defined by: 

</i|/ 2 ) = — / (/ GM )- 1 / 1 (c)/ 2 (c)d 3 «. (7) 
no J 

The spectrum of A is analytic in k = and it can be shown to 
contain a D + 2-fold degeneracy at the origin, corresponding 
to local conserved quantities. In order to solve the eigenvalue 
problem associated with (6), it is worth first to attempt the 
analysis of the long-wavelength limit k — > 0: 



L9i(c) = Xi*i(c). 



(8) 



The operator L is found to be symmetric and negative 
semidefinite with respect to the scalar product (7), hence 
eigenfunctions are orthogonal and form a complete set. In par- 
ticular, a subset of them, which spans a (D + 2)-dimensional 
subspace of the Hilbert space can be found corresponding to 
the degenerate zero eigenvalue. These are the collision invari- 
ants / GM X°, with X° denoting a set of lower order Sonine (or 
associated Laguerre) polynomials: 



X° = 



1,2c, 



C 2 



(9) 



(see also Eq. A2a). A perturbative approach is followed in or- 
der to extract those eigenvalues, denoted hereafter by Whydro, 
which reduce to zero in the long wavelength limit, from the 
full spectrum of A. The yet unknown eigenfunctions and 
eigenvalues are expanded in powers of the wavevector k: 

l*a) = \*£ ) ) + k\*U) + k 2 \*W)+..., 

u a - 4 0) + kcjW + k 2 u® + . . . , (10) 

where l^i ') denotes a linear combination of the eigenfunc- 
tions of the unperturbed system. The result of this standard 
procedure is a polynomial expression for the set {ui a } of hy- 
drodynamic modes, up to second order: 

L 

u 2 = -icofc-fc 2 (^ 0) |( C:c + co)i(c x + co)|^ 0) ), 



cj 3 = -Ar(* 



2/ lT ,(0)| 



cj 4 = -k z (-$ 



2AT,(°)l 



-cJ* 



(0)\ 



u 5 = -fc 2 (*f |C X JC X \^), 



(11) 



where cq = (5fcBTb/3m) 2 is the speed of sound of an ideal 
gas. 
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5f 



_ jGM + 



+ 



_ jGM ^ 



-fOM,„ 



+ / GM 5X • x 



/ GM AX ■ x 
A/ 



TABLE I: Notation used in this manuscript. Terms have been 
grouped and abbreviated as depicted in this table. / GM and f LM de- 
note global and local maxwellian, respectively, and A/ and 5f their 
"distance" from /. The third row informs about the closure discussed 
in this manuscript, while x is a set of lower order moments of /. 



B. Linear Hydrodynamics 

We denote by [n k ,u k ,T k ] the Fourier trans- 
forms of the hydrodynamics fields, for instance: 
rik = J_° dt dre~ ut ~ lW r dn(r,t), where Sn(r,t) 
is the fluctuation at time t and point r of the local particle 
number. The equations of hydrodynamics considered in [6], 
are the linearized Navier-Stokes-Fourier (NSF) equations, 
which represent balance equations for particle number 
density, momentum and kinetic energy endowed with specific 
constitutive equations for the stress tensor and heat flux: 



Lun k (uj) = 
ujT k (uj) = 



-in k • Ufe(w), 



(12a) 



.k fdP\ . k fdP\ 

-5- n k {u)-i — — T k (w) 
n \dnJ T n \oTJ n 

-^-fc 2 u fe H - — ((+ In) k • u fe (wXl2b) 
n Q n \ 3 / 



1 T °i, t \ 
-i — — k • Ufc(cj) 



no C v 



C v n 



k 2 T k (io), (12c) 



where £ and n are respectively the bulk and shear viscosity, 
C v is the specific heat at constant volume and k is the ther- 
mal conductivity. Solving the Eqs. (12) amounts to determine 
the eigenvalues of a 5 x 5 non Hermitian matrix, which repre- 
sent the decay rates of the collective excitations. The intuition 
enlightened in the paper [6] was to put into correspondence 
the macroscopic eigenvalues with their microscopic counter- 
part, obtained from (8) by application of perturbation theory. 
This identification allowed to find an approximate expression 
for transport coefficients only in terms of the one-body dis- 
tribution function which turned out to be equivalent to re- 
duced expressions determined by many-body autocorrelation 
functions. These coefficients properly recover the Chapman- 
Enskog expressions from classical kinetic theory. Within the 
above construction it is found that the decay rates of hydro- 
dynamic modes in the NSF approximation are quadratic in 
the wave vector Re(cj) oc — k 2 and unbounded. The use of a 
suitable projector, on the other hand, outlined in the next sec- 
tion, allows us to find proper asymptotic s and paves the way 
to solve the eigenvalue equation (6) as well as to determine 
exact transport coefficients. 



III. THE INVARIANT MANIFOLD TECHNIQUE 

The notion of invariant manifold is a generalization of nor- 
mal solution in the Hilbert and Chapman-Enskog method. 
Given a dynamical system 



(13) 



where J( f) is the vector field which induces the motion in the 
space of distribution functions U. Given bounded and smooth 
functions x(r, t) we define the locally finite-dimensional man- 
ifold fl C [/as the set of functions /(x(r, t),c). Hence, 
we will only consider sets of distribution functions whose de- 
pendence upon the space variable r is parameterized through 
some "moments" x(r, t). As it will be discussed in Sec. IV, 
once we identify such coarse-grained fields x(r, t) with the 
hydrodynamic fields, postulating their existence corresponds 
to invoking the hypothesis of local thermodynamic equilib- 
rium. Hence, the extent of our predictions is inherently re- 
stricted to lenght scales wherein the concept of a field as en- 
semble average over a statistically significant number of parti- 
cles is still meaningful. Let us denote by T w the tangent space 
to the manifold O at the point w of the phase space, and let 
us introduce a projection operator P which, when acting on 
J(f), describes the motion of the vector field along the mani- 
fold. The dynamics is, hence, splitted into a fast motion on the 
affine subspace w + ker[P] and a slow motion, which occurs 
along the tangent space T w [4] . The set of eigenvalues Whydro 
is determined as follows: 

1 . We seek for an invariant manifold O C U such that the 
following Invariance Equation (IE) is fulfilled: 



(1 - P)AA/ = 
where Af = f - f GM (cf. also Tab. I). 



(14) 



2. After determining the nonequilibrium distribution func- 
tion from Eq. ( 1 4), we derive equations of linear hydro- 
dynamics via integration of the kinetic equation (3). By 
construction, the decay rates of the macroscopic excita- 
tions then coincide with Whydro- 

Let x = [n, u, T] denote the set of dimensionless hydro- 
dynamic fluctuations: n = (n — no) /no = (particle num- 
ber perturbation), u = u/vt = (velocity perturbation) and 
T = (T— To) /To (temperature perturbation). Further, we split 
the mean velocity u uniquely as u = u ll e n +u ± e ± , where the 
unit vector e y is parallel to k, and e ± orthonormal to e y , i.e., 
ej_ lies in the plane perpendicular to k. Due to isotropy, u ± 
alone fully represents the twice degenerated (shear) dynamics. 
By linearizing around the global equilibrium, we write the lo- 
cal Maxwellian contribution to fin (2) as / LM = / GM (l+^ ) 
where ipo takes a simple form, ipo = X° • x (linear quasi equi- 
librium manifold), where X° (c) was defined in Eq. (9). It is 
conveniently considered as four-dimensional vector using the 
four-dimensional version x = [h, u" , T, u x ], and is then given 
by (A2a). It proves convenient to introduce a vector of veloc- 
ity polynomials, £(c), which is similar to X° and defined by 
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(A2b), such that (f GM ^\f GM X°) = 8^. Hence, the fields 
x are obtained as (£(c)) jlm = x, where averages are defined 
as: 

(t(c)) f = -W ^(c)/(c)d 3 v = (f GM m\f) ■ d5) 

We introduce yet unknown fields <$X(c, k) which character- 
ize the part Sf of the distribution function. As long as de- 
viations from the local Maxwellian stay small, we seek for a 
nonequilibrium manifold which is also linear in the hydrody- 
namic fields x themselves. Therefore, we set: 

6ip = SX • x. (16) 

The "eigen"-closure (16), which formally and very gener- 
ally addresses the fact that we wish to not include other 
than hydrodynamic variables, implies a closure between mo- 
ments of the distribution function, to be worked out in de- 
tail below. By using the above form (16) for Sf = f GM S(p, 
with LSf = / GM L[(5X] • x, and the canonical abbreviations 
AX = X°(c) + £X(c, k), Eq. (6) reads: 

cjf GM AX • x = AAf = -ik • cf GM AX • x + f GM L5X • x. 

(17) 

The microscopic projected dynamics is obtained from (14) by 
introducing the thermodynamic projection operator, defined 
in [4], which, when acting upon </(/) = AAf, gives: 

PAAf = D x Af ■ J £(c)AAfd 3 v, (18) 

where D x Af = dAf /dx and the quantity inside the inte- 
gral in (18) represents the time evolution equations for the 
moments x. These are readily obtained by integration of the 
weighted (6) as 

«<€(c))y = -ik • (€(C)C)/ + (|(C)) tsf . (19) 

As shown in Tab. I , D % Af = / GM AX holds, whereas (19) is 
linear in x and can be written as ojx = M • x. Hence, Eq. (18) 
attains the form: 

PAAf = / GM AX ■ M x. (20) 

In the derivation of (20), one needs to take into account that 
(£(c))sf = (as the fields x are defined through the local 
Maxwellian part of the distribution function only) and that 
(£( c )) tsf = 0' The dependence of the matrix elements of 
M upon moments of Sf is explicitly given in Tab. II. Com- 
bining (17) and (20), and requiring that the result holds for any 
x (invariance condition), we obtain a closed, singular integral 
equation (invariance equation) for complex-valued SX, 

AX • M = -ik • c AX + LSX. (21) 

Notice that <5X = AX — X° vanishes for k = 0, which im- 
plies that the invariant manifold flk^o m that limit is given by 
the set of local Maxwellians / LM . The implicit equation (21) 
for SX (or AX, as X° is known) is identical with the eigen- 
closure (16), and is our main and practically useful result. The 
Bhatnagar-Gross-Krook (BGK) collision model treated in [8] 
is recovered for L(5X) = —SX. 



A. Solving the Invariance Equation 

The invariance equation (21) as well as some symmetry re- 
lations for the components SX^ of the nonequilibrium distri- 
bution functions (worked out in Appendix B for the interested 
reader) are exact. Solutions to this equation can be obtained 
in simple cases. Considering the BGK kinetic equation, for 
instance, the IE could recently be solved numerically and the 
spectrum of hydrodynamic modes at arbitrary wavelength has 
been successfully determined [8]. In the present case, our 
strategy to solve (21) is to confine ourselves with a special 
kind of interaction potential (Maxwell molecules) and is based 
upon the results obtained by Chang-Uhlenbeck [9]. They pro- 
vided an analytical solution to the eigenvalue problem for the 
BE with the Maxwell molecules collision operator (i.e.: gas 
molecules interacting via a potential V oc r~ 4 , see also [10]). 
Their analysis showed that due to the isotropy of the operator 
L (i.e. it commutes with rotation operators in velocity space), 
it admits the following set of eigenfunctions ^> r ,i(c), 

r\(l + ±)0F , (r) „ 

* r ' l = \l (i + r + ±y, cPli z) 5 4( c )' w 

where Pi and S l+ i denote, respectively, Legendre and Sonine 
polynomials, c = |c| and z = c ■ e n /c (see also Appendix A). 
These eigenfuctions are orthonormal with respect to the scalar 
product (7), with corresponding eigenvalues: 

X r .i =2tt J sin(i?)F(i?)T r ,i(i?)di?, (23) 

where the explicit expressions for the F(z?) and T r i({f) are 
needed to numerically solve (21) and hence delegated to Ap- 
pendix A. Whereas the construction outlined in Sec. II de- 
duces the eigenvalues of the perturbed system (8), which van- 
ish in the k — > limit, just from the knowledge of ker[L] (i.e., 
the "ground states" of the unperturbed system), here we at- 
tempt a different route. We introduce, first, a decomposition 
of the microscopic particle velocity, where its components can 
be expressed through the absolute value of velocity, c, and the 
cosine of the angle between velocity and wave vector, denoted 
as z, see (Al). Next, we expand our functions [X(°\ c5X] in 
terms of the orthonormal basis ty r i = 5 , rj i(c, z): 

N 

Xf\c,z) = J2 a{ ° )(r ' l) *r,l(c,z), (24a) 

r,l 
N 

5X^(k,c,z) = J2 a £' l) (k)*r,l{ c , z )- ( 24b ) 

r,l 

The equilibrium coefficients are known, and can be deter- 
mined, by taking advantage of the orthogonality of the eigen- 
functions, as: 

a (o)W) =7r -| J e -- 2 ^ r l ( CjZ )X^(c,z)d 3 c. (25) 

Inserting (24) into the IE (21), we obtain the following non- 
linear set of algebraic equations for the unknown coefficients 
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(4 r ' /) + 4 0)(r '°)^= ( 26 ) 

N N 



with: 



r,J 



(r,l,r',V) 
(r,l,r',l') 



(f^ r , tl ,\L\f^ rtl ) 

(/ GM * r ,|c|/ GM * r ,,,). 



(27a) 
(27b) 



For any order of expansion, the solutions of (26) character- 
ize an invariant manifold in the phase space. The matrix ele- 
ments C^i tT i t ii\ can be easily evaluated in few kinetic models, 
as for the BGK collision operator, hard spheres and Maxwell 
molecules. In particular, the latter case is recovered by setting: 



/-•Maxw 

J ~(r,l,r'A') 



(28) 



Furthermore, the simplest case is BGK where all nonvanish- 
ing eigenvalues attain the constant value: Abgk = — 1. 

The calculation of the coefficients central in our 

derivation. Through these coefficients, the invariant manifold 
il C U is fully characterized: that is, the distribution func- 
tion is determined and the corresponding matrix M of linear 
hydrodynamics is made accessible. Generalized transport co- 
efficients such as viscosity and diffusion coefficients, defined 
in the Table II, can be expressed in terms of these coefficients 
and they enter the definition of the stress tensor and heat flux 
(explicit expressions provided in Appendix A). In the regime 
of large Knudsen numbers the coefficients a^ r '^ may be fur- 
ther used to, e.g., directly calculate phoretic accelerations onto 
moving and rotating convex particles [11]. 



B. Hydrodynamic modes and transport coefficients 

With M at hand, the hydrodynamic modes can finally be 
obtained from (20). The damping rates of the fluctuations 
(given by the real part of the hydrodynamic modes) are ob- 
tained by truncating the series (24) at the 4th order, and rep- 
resented in Fig. 2. The first important finding is that, for any 
finite order of expansion, the modes extend smoothly over all 
the wavevector domain and, for large k, they attain an asymp- 
totic value. This reflects the fact that, below a certain length- 
scale (more specifically, for lengths less than the mean free 
path), we reach the free-streaming limit, i.e., the regime in 
which the collisions cease to occur and particles move along 
straight lines. Hence, when reducing further the length scale, 
we may not expect an increase of the damping rate without 
the "thermalizing" effect of collisions. These physical argu- 
ments were already supported by the study of the BGK kinetic 
equation [8], wherein the hydrodynamic modes, in the limit of 
small wavelengths, reach all the same value equivalent to the 
constant eigenvalue of the BGK collision operator. A further 
indication of the role played by the spectrum of L for large k 
is provided by the observation that, when taking into account 



II 


°\ 


|| 

^3 


<7 4 


/ \ II X V \ 


/ \ II X V \ 


/ \ II X V \ 
{A"0X3) 


(C|| c^OYi) 


— K t> 


ikA 




IKJJ 


real, © 


imag,© 


real,© 


imag,© 
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Qi 


II 

q. 


II 

93 






(j ll 6X 2 ) 


<7"<5*3> 


<(c 2 - |)c^y 4 ) 


ikX 


-k 2 Z 


ikY 


-k 2 U 


imag,e 


real,© 


imag,© 


real,© 



TABLE II: Symmetry adapted components of (nonequilibrium) 
stress tensor a and heat flux q, both introduced in (A5). Row 2: 
Microscopic expression of these components (averaging with the 
global Maxwellian). Short-hand notation used: A" = c 2 — ^- and 
7" = (c — !)cy. Row 3: Expression of the components in terms 
of (as we show, real- valued) functions A-Z (see text). Row 4: Par- 
ity with respect to z - symmetric (ffi) or antisymmetric (©) - of the 
part of the corresponding SX entering the averaging in row 2, and 
whether this part is imaginary or real-valued (see Fig. 1). Row 3 is 
an immediate consequence of row 4. 



all the set of the eigenvalues of L which are unbounded below, 
also the hydrodynamic modes grow unboundedly. 

Generalized transport coefficients are obtained by the non- 
trivial eigenvalues of — fc 2 Re(M): A 2 = —A (elongation vis- 
cosity), A3 = — |y (thermal diffusivity) and A4 = —D (shear 
viscosity). After some algebra it is possible to recast the ex- 
pression for the higher order moments in terms of time cor- 
relation functions, in order to show the connection with the 
familiar Green-Kubo expressions. To this aim, we first write 
the nonequilibrium distribution function at time r as: 



$/(k,c,r) 



At 



<*/(M,0). 



(29) 



Then, due to (24), by integrating both sides of (21), we find: 

r,l 

+ ^ a (c) W)(/ GM^| Ae Ar| / GM^ i)) {3Q) 



where al and al are, respectively, real and 
imaginary-valued coefficients and (/ GM ^ M |Ae Ar \f GM ^ r ,i) = 
— ik(f GM S lfl |c||e AT \f GM ^ r ,i) because lower order moments of 
the collision operator identically vanish. Next, using the op- 
erator identity: 



Ae AT = -A[ / e JW dt\A 



.At, 



(31) 



we find for the real part of the M matrix in (30), for an arbi- 
trary time t = 0: 

N poo 
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FIG. 1: (Color online) All contributions 5Xi_4(c, k) vs. c (horizontal, c = |c|) and z G [—1, 1] (vertical axis, z is the cosine of the angle between k and 
peculiar velocity c) to the nonequilibrium distribution function Sf = / GM 5X M x M (16) at k = 1, obtained with the fourth order expansion, N = 4. Shown 
here are both their real (top) and imaginary parts (bottom row). 




FIG. 2: Hydrodynamic modes u> of the Boltzmann kinetic equation with 
Maxwell molecules collision operator as a function of wave number k. Shown 
are two complex conjugated acoustic modes oj ac , twice degenerated shear 
mode w s h and a thermal diffusion mode oJdiff- 



where $! r> i = A^;. Equation (32) extends to arbitrary wave 
vector the Green-Kubo relations for transport coefficients. 
These relations hold in the hydrodynamic regime, when the 
system, as a result of many collisions, has reached local equi- 
librium. The opposite regime (k ^> 1) is represented by a 
simple gas of noninteracting point particles. Importantly, as 
it is evident from Fig. 3, and as already noticed in [12, 13], 
the transport coefficients vanish in the limit of small wave- 

(r)(r I) 

lengths. This is due to the fact that the coefficients al 

(c)(r I) 

and a v , solving (26), vanish in that limit. This vanish- 
ing character of transport coefficients (and, hence, of the heat 
flux and the stress tensor as is evident from Tab. II) for large 



FIG. 3: (Color online) Moments A-Z of the distribution function (see Tab. II 
and Eq. A6) vs. wave number k obtained with the solution of (21). Non- 
triangles (black symbols): Moments entering only the longitudinal compo- 
nent of hydrodynamic equations. Triangles (blue symbols): Moments enter- 
ing the transverse component of hydrodynamic equations. 



k, corresponds to Eulerian (inviscid) hydrodynamics. We are 
led, then, to similar conclusions to those traced when we dis- 
cussed, in Sec. Ill, the k — > limit of the invariance equa- 
tion (14): in the free streaming regime, the local equilibrium 
manifold (local Maxwellian) becomes an invariant manifold. 
Let us recall that the Maxwellian distribution constitutes the 
zero point of the collision integral, in the sense that, in lo- 
cal equilibrium, the net flux of molecules entering and leav- 
ing an infinitesimal volume in space, due to the scattering 
processes, is zero. What we observe here is that, at a suffi- 



7 



ciently short length-scale, the distribution function reduces to 
a Maxwellian, since the contribution from the scattering event, 
again, vanishes: but now this is because collisions ceased to 
occur. 





108,0 01 



iog 10 uf 



FIG. 4: (Color online) (a) Damping spectrum, i.e., the negative imaginary 
part of k divided by frequency ui vs. the negative logarithm of lo. Results ob- 
tained in this work (by solving Eq. 21, and subsequently Eq. 20 for w(k) with 
complex- valued k and real- valued to) are compared with previous approaches 
including Navier Stokes (NS), regularized 13 moment (Regl3) [14], Grad 13 
moment (Gradl3), and experimental data presented in [15]. (b) Phase spec- 
trum, i.e, real part of k times velocity of sound co and divided by u> vs. the 
negative logarithm of a;. Again, we compare with reference results. 



IV. FINITE WAVELENGTHS HYDRODYNAMICS 

We raised the issue of the validity of the notion of invari- 
ant manifold for small wavelengths. We were able to show 
that hydrodynamic modes and the generalized transport coef- 
ficients extend smoothly over all the fc-domain (there is no 
occurrence of any critical point as in the case of a Grad ki- 
netic system, studied in [16]); hence our approach, here and 
in [8], tends to predict that the notion of invariant manifold 
holds also for short length scales. This would be in agreement 
with the celebrated papers by Alder etal. [ 1 2, 1 7] who consid- 
ered a hard spheres gas and showed that hydrodynamic laws 
remain valid down to times comparable with the time between 
collisions, t co \\, and that the fc-dependent zero frequency trans- 
port coefficients decay until they vanish at short length scales. 
It would be significant, therefore, to investigate the features 
of our model at finite frequencies and wavelegths and verify 
whether the procedure of truncation we introduced in (24) in- 
troduces a length scale below which our coarse-grained de- 
scription breaks down. In Fig. 4 a comparison is shown about 
inverse phase velocity and damping for acoustic waves be- 
tween our results, former approaches [14, 18] and experimen- 
tal data performed by Meyer and Sessler [15]. As it is seen, 
our results are very close to the predictions of the regularized 
13 (Reg 13) moments method [14] and closer to experimental 
data than Reg 13 concerning the phase spectrum. Our theory 
is capable to predict that the phase speed remains finite also at 
high frequencies, a feature which is not possessed by any hy- 
drodynamics derived from the CE expansion. A further clue 



about the features of our predictions in the regime of finite 
frequencies and wavevectors can be achieved by a closer in- 
spection upon the spectrum of density fluctuations. To this 
aim, we introduce the Laplace transform of the hydrodynamic 
fields x k (z) = / °° e~ zt Xk(t)dt and write the equation of lin- 
ear hydrodynamics, analogously to (20), as: 



x k (z) = (zI-M)- 1 x k (t = 0). 
By inverting the Laplace transform one obtains: 



x k (t) 



1 



01 - M) 



dz x k (i = 0). 



(33) 



(34) 



In order to proceed further and calculate the intermediate scat- 
tering functions C x . x (k, t) = (x k (i)x_ k (0)) it is needed, then, 
to define the averages which are employed in the calculation 
of correlation functions. These are, in fact, no longer ensem- 
ble averages, as in (15), but, due to (34), are averages over 
initial conditions, weighted by the probability density of ther- 
modynamic fluctuation theory [19, 20]. Finally, the power 
spectrum of C XjX is given by its Fourier transform: 



C*,x(M)e- 



(35) 



It is worth focusing upon the spectrum of density fluctua- 
tions, Sn t n, as, being it related to the scattering cross-section, 
it is a quantity which is experimentally accessible. The calcu- 
lation of Sn.n proceeds along the lines indicated above. It just 
suffices to notice how the solution for n k (z) involves terms 
proportional to the initial values of n k ,u k , T k , but, following 
standard recipes [19], only the term proportional to n k (£ = 0) 
needs to be retained in the calculation. By considering just the 
lower order terms in k, one obtains: 



n k (t) 



- X k 2 t 



-rk-'t 



10 



cos(cofci) 



n*(0). (36) 



The first term in (36) represents a fluctuation which de- 
cays according to a purely diffusive process, with a lifetime 
proportional to Dt, whereas the second term represents a 
fluctuation propagating through the fluid at the (dimension- 
less) speed of sound co = \/5/3 and decaying with a life- 
time given by V. The coefficient Dt generalizes the standard 
thermal conductivity, while T generalizes the combined ef- 
fect of both thermal conductivity and longitudinal kinetic vis- 
cosity. In the limit of small k, and following standard text 
books [21], their expression is given by Dt = ^(X — Y) and 



r = \X + jgY). Unlike standard treatments of hy- 

drodynamic fluctuations, the generalized transport coefficient 
X enters the expression of the coefficients Dt and V, even 
though its contribution, as it is evident from Fig. 3 is fairly 
small. The (approximate) intermediate correlation function is 
then obtained by averaging: 



Cn, fi (k,t) = 4,o + (n k (0)n- k (0)) x 



2 -D T k*\t\ + 3 -r^| t | cos(coH) 
5 10 



(37) 
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and the dynamical structure factor, hence, attains the follow- 
ing form: 

Sn,n(k,uj) = 5(cj)5 k , + <n fc (0)n_ fc (0)> x (38) 
"2 2D T k 2 3 2Tk 2 

lu 2 + {D T k 2 ) 2 + 10 (uj±c k) 2 + (Tk 2 ) 2 _ ' 




FIG. 5: (a) Dynamic structure factor S(k, u>) vs. w for a small k = 0.4 and 
(b) large k = 100. ui a = cgk denotes the hydrodnamics predicted sound 
mode of the spectrum, and the widths are related to the moments A— Z (see 
Fig. 3). For small k, these are given by Dt = \(X-Y)m&T = -{\A+ 
■|X + j^y), where A is the generalized longitudinal kinetic viscosity, Y the 
generalized thermal diffusion coefficient and X is a cross-coupling transport 
coefficient, relating heat flux to density gradients, (c) Width Dj<k 2 of the 
Rayleigh peak vs. k (double-logarithmic). At small k, D^k 1 oc fc 2 as 
all moments A-Z, except X, reach a finite value in this limit. The inflection 
point at k = k* (N) 2> 1 (shown to be increasing with the order of expansion 
N) denotes the onset of departure from the ideal Maxwellian behavior, where 
the width of the peak starts to behave sublinearly in k, and is used to quantify 
the range of validity for results obtained at finite order. 

Representative plots of S(k,u>) are shown in Figs. 5a-b. 
For small k (hydrodynamic limit), the spectrum we obtain 
recovers the usual results of neutron (or light) scattering ex- 
periments and consists of three Lorentzian peaks. The one 
centered in lu = is the Rayleigh peak, which corresponds 
to the diffusive thermal mode. The two side peaks centered 
in u> ± Crjfc are the Brillouin peaks, and represent the two 
propagating sound waves. By increasing the wave-vector, the 
structure of (38) is unchanged except that the generalized co- 
efficients Dt and T need to be replaced by more complicate 
expressions, not given here. The net effect observed is that 
sound waves get strongly damped and vanish, whereas the 
central Rayleigh peak decreases and broadens. Density fluctu- 
ations are, therefore, driven only by a diffusive thermal mode 



for large enough k. A deeper look about the behavior of the 
width at half maximum of the central Rayleigh peak with in- 
creasing wavevectors allows us to bridge the gap between the 
hydrodynamic continuum-like description and the free parti- 
cle limit. The hydrodynamic regime is featured by a width 
increasing with the square wavevector, oc k 2 . On the con- 
trary, in the free particle limit, the calculation of the dynamical 
structure factor S n:n (k,uj) reduces to the Fourier transform 
of the self part of the van Hove function G s (r, t)[19], which, 
upon writing c = r/t, is given by the Maxwellian distribution: 
G s (r,t) = Tr- 3 / 2 VT 3 t- 3 exp(-r 2 /t 2 ). Hence, the width of 
the peak is expected to grow up linearly in k, for large k. Our 
results, see Fig. 5c, predict a width which is truly quadratic 
for small enough k, reach the regime of linear behavior and 
terminate, for some large k, with a sub-linear dependence on 
k. The onset of the terminal regime at k — k*(N) marks the 
range of validity which can be accessed at a given finite order 
of expansion, N. Increasing N thus does not alter the over- 
all picture we obtained at a moderate order of expansion, and 
more generally, results obtained with N + 1 will not change 
those obtained with N below k*(N), cf. Fig. 5c. 



V. CONCLUSIONS 

The main result of our paper is the characterization of the 
nonequilibrium distribution function, through the method of 
invariant manifolds, and the calculation of its moments (the 
functions A-Z), which constitute the building blocks of the 
generalized hydrodynamic equations. As we had previously 
shown in [22], the latter equations are stable and hyperbolic 
for arbitrary wavevectors. Moreover, we have proposed and 
applied a route to solve the eigenvalue problem associated 
with the BE (6), by calculating the hydrodynamic modes, 
which we may regard either as decay rates of hydrodynamic 
fluctuations as well as generalized eigenfrequencies of the BE 
(3). The generalized transport coefficients have been numer- 
ically determined and settled into expressions recovering the 
Green-Kubo formulas. Finally, also by comparing with avail- 
able experimental data and previous approaches, we discussed 
the range of validity of our approach, which turned out to 
be capable of extending the hydrodynamic scenario to length 
scales below the mean free path. This offers new perspectives 
towards a deeper comprehension of the transition between a 
"mesoscopic" particle-like description of matter and the "con- 
tinuum" macroscopic one. 
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(Cy — -C 2 )5X„d 3 C OC 1 — 8n,4,, 



c«c±SX n d c oc S n .4, 



C (X 1 — 5 n 



I C]l (c 2 -^)SX n d 
J c±(c 2 — ^)8X n d 3 c oc S. 



(39) 



APPENDIX A: OBTAINING THE IN VARIANCE EQUATION 

In order to calculate the averages occurring in Sec. Ill, like 
(/i(c)) f, we switch to spherical coordinates. For each (at 
present arbitrary) wave vector k = ke n , we choose the coordi- 
nate system in such a way that its (vertical) z-direction aligns 
with e y and that its .x-direction aligns with e^. The velocity 
vector we had been decomposed earlier as u = u ll e n + u ± e ± . 
We can then express c, over which we are going to perform 
all integrals, in terms of its norm c, a vertical variable z and 
plane vector (azimuthal angle ■ = cos <fi; the plane 
contains e ± ) for the present purpose as: 



c/c= \/T^ 



(Al) 



as shown in Fig. 6. 

The local Maxwellian, linearized around global equilib- 
rium, takes the form: / LM // GM = 1 + (p = 1 + X° • x, 
where the four-dimensional X°, and the related vector £, em- 
ploying four-dimensional x = [h, it" , T, u 1 ], are given by 



X°(c) = (l,2 C|l) (c 2 --),2 Ci ), 
l ) c n ,^(c a -|), Ci ] 



Here, we introduced, for later use, the abbreviations 
en = c-e,|, c ± = c-e ± , = c-e^ = — — — , (A3) 
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e h k/k 




e i x 611 



FIG. 6: Schematic drawing introducing a orthonormal frame ep, e_i_, and 
ej_ X e [| which is defined by the wavevector k || ey and the heat flux q (not 
shown), which lies in the ey-e^-plane. Shown is the velocity vector c (Al) 
relative to this frame (characterized by length c, coordinate z, and angle <j>) 
and its various components. The integration over d 3 c = c 2 dcdzd<f> is done 
in spherical coordinates with respect to the local orthonormal basis. 



such that ik ■ c = ifccy. We can then rewrite (Al) as c = 
c^e^ + c n e M with c y = cz and = cy/1 — z 2 . The latter 
two components, contrasted by c ± (and e^), do not depend 
on the azimuthal angle. We further introduced yet unknown 
fields (5X(c, k) which characterize the nonequilibrium part of 
the distribution function, Sip = 5f/f GM . By analogy with the 
structure of the local Maxwellian, those are linear in terms of 
the hydrodynamic fields x themselves, 

Sip = SX ■ x = JXin + SX 2 u 11 + 5X 3 f + ftti 1 . (A4) 

The functions 6X1^,3, which are associated to the lon- 
gitudinal fields, inherit the full rotational symmetry of 
the corresponding Maxwellian components, 6X12,3 = 
SX\ t 2g(c,z), whereas 6X4 factorizes as SX&(c,z,4>) = 
2SY±{c, z) J2m=i V™ cosnufr. In this context it is an impor- 
tant technical aspect of our derivation to work with a suitable 
orthogonal set of basis functions (irreducible tensors, cf. [26], 
for models beyond the Maxwell gas) to represent 5f uniquely. 
The matrix M in (20) contains the non-hydrodynamic fields, 
the heat flux q = (c(c 2 — and the stress tensor cr = 
( cc 1 ) f, where I 1 denotes the symmetric traceless part of a 



tensor s [16, 23, 24], s 1 = ±(s + s T ) - §tr(s)I. Using 
(16) and the above mentioned angular dependence of the <5X 
functions (the only term in 6X4 playing a role in our calcu- 
lations is the first order term cos(j>, with yi = 1, see [27]), 
constraints, such as the required decoupling between longi- 
tudinal and transversal dynamics of the hydrodynamic fields, 
are automatically dealt with correctly when performing inte- 
grals over <j>. More explicitely [27], the stress tensor and heat 
flux uniquely decompose as follows 



cr 

q 



II 3 ■ 

2 e " e " 

g"e,, +q ± e A 



a 2 e,,e 



(A5a) 
(A5b) 



with the moments cr 11 = (a\,a\,a\) ■ (n,u ll ,T) and = 
(T4U" 1 -, and similarly for q (see Row 2 of Tab. II). The pref- 



actors arise from the identities 

■i. We note in passing that, while the 



= I and 



stress tensor has, in general, three different eigenvalues, in 
the present symmetry adapted coordinate system it exhibits 
a vanishing first normal stress difference. Since the integral 
kernels of all moments in (A5) do not depend on the az- 
imuthal angle, these are actually two-dimensional integrals 
over c £ [0, 00] and z E [—1,1], each weighted by a com- 
ponent of 27rc 2 / GM i5X. Stress tensor and heat flux can yet be 
written in an alternative form which is defined by Row 3 of 
Tab. II. As we will prove below, due to fundamental symme- 
try considerations, the hereby introduced generalized trans- 
port coefficients A-Z are real-valued. They can be expressed 
in terms of the moments of the distribution function, i.e., ex- 
pansion coefficients a( r >'), as follows: 



A = 

C = 

Y = 

D = 

U = 



■ (0,2) 

J ' 2 ) 



(0,2) 

V3fc 2 ' 



2k 

. N 



, z = 



2k 
V5< 



,(1,1) 



2fc 2 



V* „( r >') / fGM - „ lf GM iT( \ 

r,l 

1 N K 

-^E a i r '^/ GM ( c2 -|)^i/ GM ^)- ( A6 > 



r,l 



We proceed by using these functions A-Z to split M into 
parts asM = Rc(M) - i Im(M), 



M = fc 



/ \ 
A 



\X %Y 
D 



— ik 



/o 1 o\ 

B C 
z 




(A7) 



B = 



k 2 C, and 



2 k 2 B, C — 2 
The checkerboard structure of the matrix 



with abbreviations 
Z = \{\-k 2 Z). 
M (A7) is particularly useful for studying properties of the 
hydrodynamic equations (20), such as hyperbolicity and sta- 
bility (see [16] and below), once the functions A-Z are ex- 
plicitly evaluated. We remind the reader that we use orthogo- 
nal basis functions (irreducible moments, cf. Tab. II) to solve 
(21). In order to show how the above functions enter the def- 
inition of the M matrix, we first notice that its elements are 
- a priori - complex valued. We wish, then, to make use of 
the fact that all integrals over z vanish for odd integrands. To 
this end we introduce abbreviations (0) for a real-valued 
quantity which is even (odd) with respect to the transforma- 
tion z — > —z. One notices X° = (©, 0, ©, ©), and we recall 
that A-Z are integrals over either even or odd functions in z, 
times a component of <5X (see Tab. II). 
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Let us prove the consistency of the specified symmetry of 
M and the invariance condition: Start by assuming A— Z to be 
real-valued functions. Then M^ u = © if fj, + v is even, and 



otherwise. This implies SXi 



ie, 5X 2 



© + iQ, SX 3 = Q + iQ, and SX 4 = Q + iQ, i.e., different 
symmetry properties for real and imaginary parts. With these 
"symmetry" expressions for X°, <5X, and M at hand, and by 
noticing that symmetry properties for <5X take over to L(5X) 
because the ip r j are (i) symmetric (antisymmetric) in z for 
even (odd) I and (ii) eigenfunctions of L, we can insert into the 
right hand side of the equation, L(6X) = (X°+SX.)-(M+iQ 
I), which is identical with the invariance equation (21). There 
are only two cases to consider, because M has a checkerboard 
structure, i.e., only two types of columns: Columns /x = 1 
and /i = 3: 8X^ = © + iQ because Mis^ = 0; Columns 
fx £ {2, 4}: SX^ = Q + iQ if M^i_3 = (which is the case 
for column 4) and © + iQ if = (which is the case for 
column 2). These observations complete the proof. 



APPENDIX B: EXACT SOLUTION TO THE EIGENVALUE 
PROBLEM FOR A MAXWELL-MOLECULES COLLISION 
OPERATOR 

Given the linearized Boltzmann collision operator: 

LSf = Jl dndc 1 a(n,g)gf GM (c 1 ) x 

[^(k,c)+593(k,c 1 )-^(k,c')-^(k, ci](Bl) 

where g = \v— vi \ is the absolute value of the relative velocity 
and a(n,g) the differential collision cross section. For so- 
called Maxwell molecules the collision probability per unit 
time is independent of the relative velocity: 



2K(M + to) 
Mm 



F(<&), 



(B2) 



where to, M are the masses of the colliding particles and 
F($), with $ € [0, 7r], is given in parametric form through 
the parameter eft <G [0,ir]: 



= 7T- 2 v /cos(20)A:(siii< ? i), (B3) 
[2 3 / 2 sini9sin(2^)]-Vcos2(^ 



cos 2 ((ft) K (sin <fi) — cos(2<ft)E(suii 



:, (B4) 



with the elliptic integrals K(x) = f^ 2 (l—x 2 sin 2 y) 1 / 2 dy, 

and E(x) = Jq^ 2 (1 — x 2 sin 2 y) 1 ^ 2 dy. Since the collision 
operator is spherically symmetric in the velocity space, the 
dependence of the eigenfunctions upon the direction of c is 
expected to be spherically harmonic. Indeed, the eigenvalue 
problem admits the following solutions: 



L[tp r! l(c, Z)] = \ r ,llprj( C ' Z )' 



(B5) 



? ^±i^^ ; (,)^( C 2 ),(B6) 



(i + r+|)! 



where S^ 1 ^ 2 (x) are Sonine polynomials, and Pi(z) are Leg- 
endre polynomials which act on the azimuthal component of 
the peculiar velocity c. The Legendre and Sonine polynomials 
are each orthogonal sets, 



Pi(z)P n (z)dz = ^ - - S in , 



2n I Se^c^S^^S^^dc 
Jo 1+2 +2 



7r(; + i+r)! 



r] vrp- 



Accordingly, the ip Ti i are normalized to unity with the weight 
factor /o(c) = 7r -3 / 2 exp(— c 2 ) (as defined in Sec. II A): 

f°° 2 

/ c 2 e~ c tp r ,l(c, z)tp r > t i'(c, z) dcdz 



TT- 3 / 2 / e- c2 ^( C )</v<J<(c)d 3 C. 



(B7) 



The corresponding eigenvalues for Maxwell molecules are 
given by: 



X r j = 2tt / sin(i?)F(#)T r i(i9)d0, (B8a) 

T r m = cos^^PKcos^+sin^^PKsin^) 

-(1 + ^io), (B8b) 

The collision operator is negative semidefinite, that is, all 
eigenvalues are negative except Ao,o, Ao.i, and Ai o which 
are zero and correspond to the collision invariants. As it 
was shown in [9], there is no lower bound for the set of 
eigenvalues. Chang and Uhlenbeck's investigation [9] of the 
dispersion of sound in a Maxwell molecules gas was based 
upon writing the deviation from the global equilibrium as: 



(ifa+Sip) = J2T r , 



{r,l}=0 Ul % ; 



ty r j (c) so that the eigenvalue equa- 



tion reduces to an algebraic equation for the coefficients a r j: 

oo 

ua r j = — ik • Qr,l,r',l'CI>r',t' + ^r,lO>r,h (B9) 

{r',l'}=0 



(BIO) 



The hydrodynamic modes for the Maxwell-molecules gas 
are found by setting to zero the determinant of the above sys- 
tem of linear equations. Within this approach, from the knowl- 
edge of the spectrum of L, it is possible to solve the eigenvalue 
problem (3) for an arbitrary number of modes, just by tuning 
the number of eigenfunctions taken into account in the ansatz 
for the nonequilibrium distribution function. The peculiarity 
of the Maxwell-molecules gas lies in the fact that at any stage 
of approximation the modes recover and extend those corre- 
sponding to lower order approximations. This method pro- 
duces results which are found to be in agreement with the CE 
expansion. 



